First-principles approach to lattice-mediated magnetoelectric effects 
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We present a first-principles scheme for the computation of the magnetoelectric response of mag- 
netic insulators. The method focuses on the lattice-mediated part of the magnetic response to an 
electric field, which we argue can be expected to be the dominant contribution in materials dis- 
playing a strong magnetoelectric coupling. We apply our method to Cr2C>3, a relatively simple and 
experimentally well studied magnetoelectric compound. 
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Magnetoelectric (ME) materials are insulators that al- 
low control of their magnetic properties by means of ex- 
ternal electric fields [3,13], thus attracting great techno- 
logical interest. Current research focuses on obtaining 
compounds with a robust ME behavior at ambient con- 
ditions. This is proving a major challenge, as progress 
is hampered by one fundamental difficulty: the scarcity 
of ferromagnetic insulators (not to mention ferromagnetic 
and ferroelectric multiferroics 3] ) with a high Curie tem- 
perature. An additional problem pertains to the magni- 
tude of the effect: the ME response is usually very small, 
reflecting the weakness of the spin-orbit interactions that 
are typically responsible for the coupling. 

Quantum calculations based on efficient schemes like 
Density Functional Theory (DFT) have proved very use- 
ful in studies of magnetic and ferroelectric materials, and 
are expected to facilitate progress on magnetoelectrics. 
Indeed, there is a growing number of DFT works tack- 
ling the search for new compounds Q and even propos- 
ing new coupling mechanisms [B|, However, we still 
lack a first-principles scheme to compute the ME cou- 
pling coefficients, something that is critical to aid the 
experimental work. In this Letter we introduce one such 
ab initio methodology and demonstrate its utility with 
an application to Cr203. 

Lattice-mediated ME response.- Computing the full 
ME response of a material would require quantum rela- 
tivistic simulations with simultaneously applied electric 
and magnetic fields. While possible in principle, such cal- 
culations constitute a great challenge even for the case of 
static fields, and it is convenient to look for simplifica- 
tions of the problem. 

In the following arguments we will consider an ideal- 
ized one-dimensional (ID) crystal displaying a linear ME 
effect, the generalization being straightforward. At zero 
magnetic field, the magnetization induced by the appli- 
cation of an electric field £ is given by: 



a 2 < ^-d-^m Q j which suggests that strong ME couplings 
will occur in materials displaying large dielectric and 
magnetic responses. On more physical grounds, one can 
argue that large ME effects will be associated to signifi- 
cant electronic hybridizations or orbital rearrangements 
induced by applied electric fields, as it is processes of that 
nature that may lead to a magnetic response via the spin- 
orbit coupling. It is then worth noting that (1) such a 
response to an electric field is typical of essentially all 
highly polarizable compounds used in applications and, 
most importantly, (2) such strong dielectric responses are 
never a purely electronic effect; rather, they are driven 
by the structural changes induced by the applied field. 
One can thus conclude that large ME effects will most 
likely be based on lattice-mediated mechanisms. 

Formally, the lattice-mediated contribution to the di- 
electric susceptibility is defined as xf att = x d ~ xtiec 
where xticc accounts for the purely electronic effect cor- 
responding to clamped atomic positions and lattice pa- 
rameters. The ME coupling coefficient a can also be de- 
composed in this way, and the discussion above suggests 
that ctiatt will be the leading contribution in materials 
displaying strong ME effects. We shall thus focus on its 
computation. 

Methodology.- The structural response of an insulator 
to a small electric field can be modeled in terms of the 
infra-red (IR) modes of the material, which are obtained 
from the diagonalization of the force-constant matrix at 
the r point of the Brillouin zone (BZ). (Working with 
small fields allows us to truncate all the following Taylor 
series at the lowest order possible.) Let us denote by 
the amplitude of the i-th IR mode, with i running from 1 
to -/Vtr, and by d the corresponding eigenvalue. Taking 
the Mj's and the applied field £ as independent variables, 
we write the energy of our idealized ID crystal around 
its equilibrium state as 



AM(£) = M{£) - M s = a£ + 0(£ 2 



(1) 



E({ Ui },£)= E +\Y j Q 



where a is the linear ME coefficient and we have in- 
cluded a spontaneous magnetization M s for general- 
ity. The magnitude of the ME response is limited by 
the magnetic (x m ) and dielectric (x d ) susceptibilities as 



n £[P s + Xc d icc£ + A7W{«i})] 



(2) 



where fio is the unit cell volume and we have included 
a spontaneous polarization V s for generality. Note we 
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have decomposed the linear part of the induced polar- 
ization into electronic (xtiec^) an( i lattice contributions, 
the latter being given by 

A7? latt = — , (3) 

where pf is the polarity of the i-th. IR mode, which can 
be obtained from the atomic Born effective charges and 
the mode eigenvector Q. From these expressions, the 
equilibrium value of Ui for a given £ is calculated to be 

Ui = ^£. (4) 

On the other hand, the assumption that our model crys- 
tal displays a linear ME effect implies 

1 Win 

AM{{u t },£) = a clcc £ + -Vp>„ (5) 

where the Ui 's are again assumed to be independent of £ 
and we have introduced p™ parameters that quantify the 
magnetic response to the IR modes. Finally, by virtue of 
Eq. ([J]) we can write 

1 ^IR m J 

AM(£) = a clcc £ + -]T P -^£ , (6) 

z=l 

and it immediately follows that 

Nm i NlR p m p d 

ttlatt = >> a latt,i = 77- >> V» ' ' ( 7 ) 

where mode-dependent contributions to ai a tt have been 
defined. This equation encapsulates our method for an 
ab initio computation of the ME response. Its most re- 
markable feature is that all the parameters that appear 
in it can be computed without the need of simulating the 
material under applied electric or magnetic fields, which 
brings the calculation of ME effects within the scope of 
the most widely used DFT codes. 

The above expression offers some insight into the mi- 
croscopic ingredients needed to have a strong lattice- 
mediated ME response. In essence, one would like to 
have soft IR modes (i.e., with small d) that are highly 
polarizable (i.e., with large pf) and cause a large mag- 
netic response (i.e., with large pf 1 ). Note that the pf 1 
parameters can be viewed as the magnetic analogue of 
the polarities pf associated to the dynamical charges. 
Also in analogy with pf, pf can be written as a sum 
of atomic contributions weighted by the corresponding- 
components of the IR eigenvector. It is then clear that, 
in order to have IR modes with simultaneously large pf 
and pf 1 , we need materials in which the magnetic atoms 
present large Born effective charges. While rare, this is 
the case of compounds like CaMn03 0]. 



A few additional comments are in order. (1) The pro- 
posed approach can be used to simulate ME effects of ar- 
bitrary order in £, but is restricted to couplings that are 
linear in the magnetic field. This limitation can be reme- 
died by simulating the material under applied magnetic 
fields, which is relatively easy in comparison with treat- 
ing finite electric fields in extended systems described 
with periodic boundary conditions 0. (2) The ME re- 
sponse mediated by the strain (77) can be trivially in- 
cluded. More specifically, the above formulas are correct 
for paraelectrics, for which the leading strain terms in 
Eq. (|5J) are of the form rf and rjuf, and result in ME 
responses of order higher than linear. In ferroelectrics, 
on the other hand, there exist r]Ui terms that give a lin- 
ear contribution to the response. (3) While the above 
derivation is made in terms of the eigenvectors of the 
force-constant matrix, one could imagine an analogous 
scheme using the IR eigenmodes of the dynamical ma- 
trix as structural variables. It would then be possible to 
model the dynamical ME response. (4) It is possible to 
derive the above results in a more general way, by work- 
ing with an energy E({ui\, £; {rrij}, 7i) that includes the 
localized magnetic moments rrij and the magnetic field 
H as independent variables of the system. 

Results for C^O^.- The work on magnetoelectrics 
starts with the prediction [1(3] and experimental confir- 
mation [ll|, [13] that linear ME effects occur in Cr203 
(chromia), which remains one of the simplest and best 
studied ME materials. Ci^Os is an antiferromagnetic 
(AFM) insulator with a 10-atom unit cell and the mag- 
netic structure sketched in Fig. Q] The magnetic easy 
axis lies along the rhombohedral direction c. The crystal 
has the magnetic space group R3'c', which preserves all 
the crystallographic symmetries of i?3c; thus, the com- 
pound is paraelectric. Cr2C>3 presents six IR modes: two 
polarized along the rhombohedral c-axis, corresponding 
to the A<i u irreducible representation of 3m, and four 
double-degenerate modes with E u symmetry and polar- 
ization within the afr-plane. The linear ME tensor is di- 
agonal with two independent terms: a aa — ctbb = ol x 
and a cc = an. Naturally, the lattice-mediated part of a± 
(a\\) can be decomposed into contributions from the E u 
(A2 U ) modes, which we can compute with our method. 
(In the following we drop the "latt" subscript from the 
a's to alleviate the notation.) 

For the calculations we used the LDA [l3} approxi- 
mation to DFT as implemented in the plane-wave code 
VASP [H| . We used the PAW scheme [lf| to represent 
the atomic cores. Only the nominal valence electrons 
were explicitly solved, which we checked is sufficient. Let 
us just note that all the trivial calculations involved in 
this study (e.g., for the equilibrium atomic structure, 
force-constant matrix, or induced polarizations [13]) were 
performed accurately and following well-established pro- 
cedures, and that all of them were done at the collinear 
level. To obtain the pf 1 parameters in Eq. [[7]). we com- 
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FIG. 1: Panel a: Primitive cell ol Cr203. Solid arrows repre- 
sent the AFM ground state. Dashed arrows sketch the atomic 
displacements within the ab plane associated to a typical E u 
IR mode, as well as the induced spin rotations that render a 
net magnetic moment. Panel b: Computed polarization and 
magnetization induced by the condensation of the IR modes. 
Dashed and solid lines correspond to A 2u and E u modes, re- 
spectively. Note that the polarizations and magnetizations 
associated to the E u (A2 U ) modes lie within the ab plane 
(along the c direction). Note also that the magnetization in- 
duced by the A 2u modes is essentially zero. 



puted the magnetic response upon condensation of the 
IR modes by running fully self-consistent non-collinear 
simulations including spin-orbit couplings. Interestingly, 
we found that a non-self-consistent approach, as usually 
employed for the computation of magnetic anisotropy en- 
ergies, renders qualitatively incorrect results in this case. 
Let us also stress that, given the small magnitude of the 
energy differences associated to the ME effects in Cr203, 
one has to be very careful with the choice of the pa- 
rameters controlling the accuracy of the calculations. In 
particular, we found it necessary to use a very demand- 
ing stopping criterion for the self-consistent-field calcula- 
tions (namely, energies converged down to 10 -10 eV) to 
obtain, in a computationally robust way, reliable values 
of the magnetic moments induced by the condensation of 
the IR modes. We also determined that a fc-point grid 
of at least 7x7x7 is needed for accurate BZ integrations. 
(A magnetic easy axis in the afr-plane is incorrectly pre- 
dicted if grids that are not dense enough are used.) The 
plane-wave cutoff was found to be less critical; we used 
400 eV. We employed the "LDA+U" scheme of Dudarev 
et al. [ll| for a better treatment of the 3d electrons of Cr. 
We chose U c s = 2 eV, which renders results in acceptable 
agreement with experiment for the atomic structure, IR 
phonon frequencies, electronic band gap, and magnetic 
moments [19(. At any rate, we checked the choice of U e g 
is not critical, even for the computation of ME coeffi- 



TABLE I: Parameters of Eq. (0 computed for the IR modes 
of Cr2 03. Modes are divided in two groups, A 2u and E u , 
according to their symmetry. The last line shows the results 
for the two independent a coefficients, obtained from the ad- 
dition of the corresponding mode contributions. The a's are 
given in Gaussian units (g.u.) [22T |. 





A 2u 


modes 




E u modes 




Gi (eV/A 2 ) 


10.8 


25.7 


10.4 


16.9 21.6 


32.5 


pf (\e\) 


0.39 


8.52 


0.65 


0.16 3.24 


7.14 


p? (io-V/A) 


0.02 


0.04 


0.41 


-2.70 11.32 


8.51 


on (10- 4 g.u.) 


0.00 


0.00 


0.01 


-0.01 0.62 


0.68 


(io- 4 g.u.) 


an = 


0.00 




a ± = 1.30 





cients. Finally, let us note the orbital degrees of freedom 
can be expected to be quenched in CT2O3; thus, we ne- 
glected their contribution to the magnetization. 

Tableland Fig.[T]summarize our results, which present 
the following features. (1) We obtain an much smaller 
than a±. Indeed, our calculations indicate that the mag- 
netic response associated to the A 2u modes is nearly zero, 
and provide an explanation for such an effect. We find 
that, for the E u modes, the induced in-plane magneti- 
zation occurs via a canting of the Cr spins, as sketched 
in Fig. Q] In contrast, in the case of the A 2u modes, 
no symmetry- allowed spin canting can induce a mag- 
netization along the c direction. Instead, the simula- 
tions show that the magnetization originates from a tiny 
charge transfer between the spin-up and spin-down Cr 
sublattices. Probably, the smallness of the correspond- 
ing pf 1 coefficients reflects the relatively large energy cost 
associated to such a mechanism. (2) The ME response 
a±_ is dominated by the hardest E u modes and, inter- 
estingly, such a result could have been anticipated from 
the mode eigenvectors. More precisely, the two hardest 
modes present a relatively large Cr contribution, which 
should lead to relatively large values of pf 1 , as we indeed 
find. In addition, in the hardest eigenmode the Cr and O 
sublattices move rigidly and in opposite directions, which 
must result in a large pf, exactly as found. (3) We obtain 
both positive (from three modes) and negative (from one 
mode) contributions to a±. (Given the smallness of the 
magnetic effects computed, we have not tried to identify 
the electronic underpinnings of having positive or nega- 
tive ai's.) This result suggests that, in a general case, a 
small static ME effect may be the result of cancellations 
between contributions from different IR modes. Hence, 
large static ME effects will most likely be associated to 
compounds in which a single IR mode dominates the re- 
sponse. 

To the best of our knowledge, the low temperature 
ME response of Cr203 is not totally understood, which 
reflects both the difficulties involved in ME measure- 
ments and the rich nature of the problem. The exper- 



4 



imental results at 4.2 K are quite scattered [23[: |aj_| 
ranges from 0.2xl0~ 4 to 4.7xl0~ 4 in Gaussian units 
(g.u.) and |a||| from 0.4xl(T 4 to 1.2xl(r 4 g.u. There 
are reasons to believe that the magnitude of the ME ef- 
fects was underestimated in the early experiments [24j, 
and that the largest coefficients measured [2^, [25[ are 
the most reliable ones. In particular, |cv_l| probably lies 
somewhere between 2xl0 -4 and 4xl0~ 4 g.u., which is 
remarkably close to our result. Interestingly, it is not 
clear how to explain this relatively large value of a± in 
terms of the purely electronic mechanisms typically con- 
sidered [25|, [2(1 H3] • It is thus worth noting our computed 
lattice-mediated ME response is of the same magnitude 
as the one measured. As for the parallel response, all the 
experiments render |a||| < at low temperatures, 

but none reports an essentially zero value as we obtain. 
Our results are thus compatible with the notion that ei- 
ther a purely electronic mechanism, as the electric-field- 
induced g shift proposed in Ref. I27I , or a magnetic effect 
not related to the ME coupling [26| is responsible for the 
non-zero a\\ at low temperatures. 

In summary, we have introduced an efficient method 
to compute lattice-mediated magnetoelectric responses 
ab initio. We hope our work will enable a more effective 
interaction between theory and experiment in the search 
for materials that can be used in applications. 
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